function [paramsOpt,theta2Adj,input] = runInduceStatVAR_threshold_macro(theta2,...
    VarU_x2,Cov10U_x2,Cov01U_x2,x2,macros,econAct,thresholdLevel,h0RegimeOn,hxRegimeOn)
dispOn = 1;
factors     = [macros;x2];
[nx,T]      = size(factors);
nm          = size(macros,1);
nx2         = size(x2,1);
VarU        = zeros(nx,nx,T);
Cov10U      = zeros(nx,nx,T);
Cov01U      = zeros(nx,nx,T);
for t=1:T
    VarU(nm+1:nm+nx2,nm+1:nm+nx2,t)   = VarU_x2(:,:,t);
    Cov10U(nm+1:nm+nx2,nm+1:nm+nx2,t) = Cov10U_x2(:,:,t);
    Cov01U(nm+1:nm+nx2,nm+1:nm+nx2,t) = Cov01U_x2(:,:,t);
end

[~,~,hx_1,hx_2] = unfoldTheta2_threshold_macro(theta2,nx);
if hxRegimeOn == 1
    % Regime 1
    eigValues_1         = eig(hx_1);
    modulus_1           = sqrt(real(eigValues_1).^2 + imag(eigValues_1).^2);
    if any(modulus_1>=1)
        %delta1 = 1/max(modulus_1)*0.99;
        params.delta1 = 1/max(modulus_1)*0.99;
    else
        %delta1 = 1;
        calibrateParams.delta1 = 1;
    end
    % Regime 2
    eigValues_2         = eig(hx_2);
    modulus_2           = sqrt(real(eigValues_2).^2 + imag(eigValues_2).^2);
    if any(modulus_2>=1)
        %delta2 = 1/max(modulus_2)*0.99;
        params.delta2 = 1/max(modulus_2)*0.99;
    else
        %delta2 = 1;
        calibrateParams.delta2 = 1;
    end
else
    % We have no regime-dependent loadings, so we just check stability
    % using hx_1
    eigValues_1         = eig(hx_1);
    modulus_1           = sqrt(real(eigValues_1).^2 + imag(eigValues_1).^2);
    modulus_2           = 0; 
    if any(modulus_1>=1)
        %delta1 = 1/max(modulus_1)*0.9;
        params.delta1 = 1/max(modulus_1)*0.99;
    else
        %delta1 = 1;
        calibrateParams.delta1 = 1;
    end
    %delta2 = 1;
    calibrateParams.delta2 = 1;
end
if ~exist('calibrateParams','var')
    calibrateParams = [];
end
%params.delta = min(delta1,delta2);

if any(modulus_1>=1) || any(modulus_2>=1) 
    if dispOn == 1
        disp('OPTIMIZATION TO INDUCE STATIONARITY OF THRESHOLD MODEL')
    end
    % Moment-based scaling to induce stationarity
    meanFactors           = mean(factors,2);
    T                     = size(factors,2);
    input.varXData        = (factors-repmat(meanFactors,1,T))*(factors-repmat(meanFactors,1,T))'/(T-1)...
        -mean(VarU,3);
    input.VarU            = VarU;
    input.Cov10U          = Cov10U;
    input.Cov01U          = Cov01U;
    input.factors         = factors;
    input.theta2          = theta2;
    input.selectParams    = fieldnames(params);
    input.calibrateParams = calibrateParams;
    input.econAct         = econAct;
    input.thresholdLevel  = thresholdLevel;
    input.h0RegimeOn      = h0RegimeOn;
    input.hxRegimeOn      = hxRegimeOn;
    input.numSim          = 1D5;
    input.burnin          = 1D4;
    maxIter               = 10000;
    tolFun                = 1e-6;
    options       = optimset('Display','iter','MaxIter',maxIter,'MaxFunEvals',10000,'TolFun',tolFun,'TolX',1e-6);
    paramsValues  = struct2array(params)';
    
    paramsValuesOpt = fminsearch(@induceStatVAR_threshold_macro,paramsValues,options,input);
    [Qopt,theta2Adj,paramsOpt] = induceStatVAR_threshold_macro(paramsValuesOpt,input);
    for i=1:size(paramsValuesOpt,1)
        name = input.selectParams{i};
        disp([name ' = ', num2str(paramsValuesOpt(i,1)), ' to induce a stationary threshold model. Qopt = ', num2str(Qopt)]);
    end
else
    paramsOpt.delta = 1;
    theta2Adj        = theta2;
end

end


